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ABSTRACT 

We assess the stability of Kruskal-Schwarzschild (magnetic Rayleigh- Taylor) type modes for accreted 
matter on the surface of a neutron star confined by a strong (> 10 12 G) magnetic field. Employing the 
energy principle to analyze the stability of short-wavelength ballooning modes, we find that line-tying to 
the neutron star crust stabilizes these modes until the overpressure at the top of the neutron star crust 
exceeds the magnetic pressure by a factor ~ 8(a//i), where a and h are respectively the lateral extent 
of the accretion region and the density scale height. The most unstable modes are localized within a 
density scale height above the crust. We calculate the amount of mass that can be accumulated at the 
polar cap before the onset of instability. 

Subject headings: stars: neutron — pulsars — accretion, accretion disks — MHD — instabilities 



1. INTRODUCTION 

Despite recent progress on the spreading of accreted 
^ ' matter from a disk over an unmagnetized neutron star 
OO \ (Inogamov & Sunyaev 1999; Popham & Sunyaev 2000), 
^vO ■ there has been scarcely any effort directed at understand- 
ing how the magnetically funneled accretion stream of an 
accretion-powered X-ray pulsar spreads over the neutron 
star's surface. The presence of a strong magnetic field 
raises the possibility that the accreted matter is confined 
to the polar regions of the neutron star. This funneling 
' and confinement of the accreted matter has been proposed 
(Joss & Li 1980) as the reason for the apparent stabil- 
q ■ ity of thermonuclear burning on strongly magnetized neu- 
' tron stars. The presence of strong magnetic fields, with 
r/3 , B > 10 12 G, is deduced from cyclotron features in their 
& ■ spectrum (for reviews, see White et al. 1995; Heindl et al. 
£>. \ 2000) and from regular torque reversals (Bildsten et al. 
• »-H . 1997) that suggest the magnetospheric and co-rotation 
' radii are similar. In neutron stars with weak magnetic 
H \ fields (as inferred by the absence of pulsations or cyclotron 
features in the persistent emission), unstable burning is 
observed as type I X-ray bursts (Lewin et al. 1995). The 
funneling and confinement of the accreted matter by the 
magnetic field is believed to increase the local accretion 
rate to near- or super-Eddington magnitudes, for which 
the burning is stable to temperature perturbations (see 
Bildsten 1998, for a review). 

Confinement of the accreted matter requires its over- 
pressure be balanced by the magnetic field. Magnetostatic 
equilibria of this type were discussed by Hameury et al. 
(1983) and Brown & Bildsten (1998). In these papers, 
however, the stability of such equilibria was not addressed. 
This question is considered here. 

If the magnetic field has an approximately potential 
dipole shape (as can be expected as long as the accreted 
plasma overpressure is low compared to the magnetic pres- 
sure), the force of gravity has a nonvanishing compo- 
nent perpendicular to the magnetic field. Such a config- 
uration might be expected to be susceptible to Kruskal- 
Schwarzschild ("magnetic Rayleigh- Taylor" ) interchange 
instabilities (Kruskal & Schwarzschild 1954). Indeed, it 



is easy to show that for characteristic parameters of polar 
cap plasmas the stabilizing effect of the magnetic curva- 
ture is small compared to the gravitational drive. Never- 
theless, the line-tying in the neutron star crust eliminates 
interchange modes. Instead, ballooning modes can be un- 
stable if the ratio of plasma overpressure to magnetic pres- 
sure A/3 is sufficiently large. This instability has been a 
subject of numerous investigations in the context of labo- 
ratory (e.g., Ohkawa & Kerst 1961; Coppi & Rosenbluth 
1966; Connor et al. 1979; Freidberg 1987), solar (e.g., 
Hood 1986; De Bruyne & Hood 1989; Strauss & Longcope 
1994), magnetospheric (e.g., McNutt et al. 1987; Hameiri 
et al. 1991) and astrophysical (e.g., Parker 1966; Shibata 
et al. 1989) plasmas. In the present context we find that 
this instability occurs when Af3 exceeds a/h where a is 
the polar cap radius and h is the density scale height near 
the neutron star crust. For such a large overpressure the 
equilibrium magnetic field is strongly distorted away from 
the potential field. 

The structure of this paper is as follows. We begin by 
reviewing (§ 2) the conditions at the polar cap of an ac- 
creting neutron star. We then formulate, in section 3, the 
mathematical model which serves as the framework of our 
linear stability analysis. In section 4 we discuss the en- 
ergy principle and derive a simple intuitive expression for 
the magnetohydrodynamic (MHD) potential energy, rele- 
vant for plasmas in the neutron star polar caps. We then 
use this expression to derive a criterion for stability of 
short-wavelength (in the direction transverse to the mag- 
netic field) ballooning modes and discuss the stability of 
a plasma confined by an approximately potential, dipole 
magnetic field (§ 5). We then derive an approximate an- 
alytic expression for an equilibrium in which a large over- 
pressure distorts the field away from the potential field 
(§ 6). Using this equilibrium we find the marginal stabil- 
ity criterion (§ 7), from which we compute the amount of 
mass that must be accreted in order to drive the instabil- 
ity (§ 8). We conclude with the discussion of our results 
in Section 9. 

2. THE STRUCTURE OF POLAR CAPS 
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We are concerned with the stability of accreted matter 
on the polar cap of a strongly magnetized neutron star. 
Figure 1 provides a schematic of the situation. The col- 
limatcd infalling matter approaches the neutron star at 
roughly free-fall velocity vg ~ c. At the near- or super- 
Eddington local accretion rates, some form of radiative 
shock is likely to form (Basko & Sunyaev 1976). The region 
we arc interested in is far below where the infalling ma- 
terial decelerates and becomes incorporated into the neu- 
tron star atmosphere. The continual deposition of matter 
into the upper atmosphere pushes the underlying material 
deeper, and the downward flow velocity is v w rh/p <C c s , 
where rh is the local accretion rate per unit area. The 
atmosphere is, to a high accuracy, in hydrostatic equi- 
librium. On time scales for matter to flow through the 
thickness of the atmosphere, the field lines co-move with 
the fluid (Brown & Bildstcn 1998). Under these conditions 
we are justified in describing the structure of the accretion 
region with a magnetostatic approximation. 

The polar cap is characterized by two length scales: its 
radius and the density scale height. Various estimates of 
the polar cap radius a (Lamb et al. 1973; Eisner & Lamb 
1977; Arons & Lea 1980) suggest that a ~ 10 5 cm, and we 
shall adopt this value throughout the remainder of this pa- 
per. The density scale height h is much smaller than this 
value of a (see below). With these approximations, the 
mass and radius of the neutron star enter only through the 
surface gravity g, for which we use g — 1.86 x 10 14 cm s~ 2 , 
which is the Newtonian value for a canonical neutron star 
of mass M = 1.4M© and a radius R — 10 km. 
Infalling matter 

I Intalling matter 




Photosphere _ 
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Ocean 
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Shock 
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Fig. 1. — Overview of the accretion region. 

Numerical solutions (Hameury et al. 1983; Brown & 
Bildstcn 1998) of the magnetostatic equations show that 
the greatest distortion of the magnetic field occurs near 
the boundaries where the field lines are assumed to be 
tied. In the outermost layers of the neutron star, such a 
condition occurs at the neutron star crust where the ions 
form a Coulomb lattice. It is the rigidity of the crust that 
supports the superincumbent mass of the accreted matter 
in the polar cap and anchors the magnetic field lines. The 
electrostatic coupling is parameterized by 

r _ Z 2 e 2 /47m N \ 1/3 
k B T \ 3 J 



1.1 



Z 2 \ 
A^J 



10 8 g cm 3 



1/3 



10 s K 



(1) 



where fee is the Boltzmann constant, p is the density, and 
T, tin, Ze, and A are, respectively, the temperature, the 
number density, the charge, and the atomic number of nu- 
clei. For r > 1, the ions form a liquid. For sufficiently 
large L, the ions form a crystalline lattice. The composi- 
tion of the deep ocean and crust is unknown (see Schatz 
et al. 1999); for definiteness we take the ocean and crust 
to be composed of pure iron, and assume that solidifica- 
tion occurs for L > 173, in accordance with recent Monte- 
Carlo simulations (Farouki & Hamaguchi 1993). At this 
location, the crust of the neutron star can support a shear 
stress and therefore anchor the field lines. In the absence 
of a full calculation of the crust and magnetic field, we shall 
treat the crust [similarly as in previous works (Hameury 
et al. 1983; Brown & Bildsten 1998)] as a rigid, perfectly 
conducting surface to which the field lines are anchored 
(see, however, § 9). 

As evident in equation (1), the location of the crust is 
determined by the temperature of the ocean. We assume 
that the burning of accreted hydrogen and helium is sta- 
ble (this assumption can be checked o posteriori by finding 
the depth at which matter spreads). Stablity of the H/Hc 
burning requires that the accretion rate be locally super- 
Eddington (see Bildstcn 1998); for the remainder of this 
work we shall use the local Eddington rate as our fiducial 
value. In this case, the temperature in the ocean is set by 
the nuclear flux released from burning of the accreted fuel 
to iron peak elements (about 8 MeV per nuclcon) and is 
> 5 x 10 8 K (Brown & Bildsten 1998; Schatz et al. 1999). 

For the density and temperature of interest, the elec- 
trons are degenerate, i.e., k B T < Ep — m e c 2 , where 
E F = m e c 2 [l + (37r 2 n e ) 2 / 3 (fi/m e c) 2 ] 1 / 2 is the electron 
Fermi energy, and m e and n e are the electron mass and 
density, respectively. The principal contribution to pres- 
sure p is the electron pressure p e and therefore to a good 
approximation p w p e w n e E-p /A oc p 1 where 7 = 4/3 and 
the proportionality constant is only a function of the ion 
charge-to-mass ratio Z/A. For a plasma of uniform compo- 
sition this constant is independent of both time and space. 
We shall use this approximation in the next section. 

Using this approximate equation of state, we write the 
density scale height h — jp/pg at the top of the crust in 
terms of T, 



h w 7 



9tt 
256 



1/3 



i N A k B T r 



Z 2 / 3 A 
= 6.90 x 10 3 cm 
26 \ 2/3 /56 
~Z) \A 



5 x 10 8 K 



L 

173 



(2) 



where 7 = 4/3, N A = 6.02 x 10 23 g" 
stant and «f = e 2 /he. 

3. MODEL 



is Avogadro's con- 



We model the polar cap as the region of the neu- 
tron star atmosphere and ocean that is bounded by the 
same magnetic flux surface within which accretion occurs. 
This region is characterized by a somewhat higher pres- 
sure than in the surrounding matter, due to the weight of 
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the accreted material. The plasma in the polar cap and 
in the surrounding region will be treated as a perfectly 
conducting fluid of uniform composition, bounded by a 
perfectly conducting, rigid surface (the top of the crust) 
and immersed in a magnetic field. The magnetic field lines 
are assumed tied to the crust; the intersection point of the 
field line with the top of the crust will be referred to as 
the footpoint. Since the density scale height of the plasma 
in the polar cap is much smaller than the stellar radius 
we shall consider the field lines as semi-infinite and having 
only one footpoint. In the equilibrium, the plasma and the 
immersing magnetic field are assumed to be axisymmetric 
with the axis of symmetry coinciding with the magnetic 
axis; the magnetic field is assumed to be untwisted. 

We describe the plasma by the single-fluid MHD theory. 
In the presence of the gravitational potential U and mag- 
netic field B the momentum balance equation for a plasma 
of density p and pressure p, moving with velocity v, is 



dv 1 

p— = -J x B - Vp - pVU. 

at c 



(3) 



While Ohm's law is modified by gravity, the induction 
equation is unaffected (cf. Bernstein et al. (1958)), so that 



~dt 



= V x v x B. 



(4) 



The above equations are supplemented by Ampere's law, 



47T 



V x B, 



the continuity equation, 

dp 
dt 



V • (pv) = 0, 



(5) 



(6) 



and an equation of state which we assume to be adiabatic: 



dt 



VP 



0. 



(7) 



As mentioned in §1, the primary contribution to the 
plasma pressure in the region of interest is due to de- 
generate, relativistic electrons (7 = 4/3). Therefore, in 
region of uniform composition, pp 1 is, to a high ac- 
curacy, constant in space. For the assumptions in § 2, 
e = dXnpp-^ydlnp w 3/8p io „/;Pc - (3/2Z){k B T/ E F ) w 
4 x 1CT 4 . Thus 

pp^ 1 = const., (8) 

to a very good approximation. 

The force balance in a static equilibrium is given by 



1 



-JxB=Vp + pVU, 



(9) 



from which it follows that 



Vp + pVU = pVF, (10) 

where the function F is a function constant on the field 
lines 

B-VF = 0. (11) 

By setting the gravitational potential to zero at the top of 
the crust, F is determined by the sound speed c s = \fjpjp 
there, F = C y( 7 -l). 



An untwisted, axisymmetric equilibrium magnetic field 
can be represented in the form 



B = Vip x V9 



(12) 



where tp is the flux function and 9 is the toroidal angle. 
It is straightforward to show with the aid of Ampere law 
that in such a geometry the current is perpendicular to 
the magnetic field: J • B = 0. Eq. (11) implies that F is a 
function of the flux function only: F — F(ip). 

Upon introducing cylindrical coordinates (r, 0, z) with 
the ^-direction along the axis of symmetry, the normal 
(i.e., parallel to VVO component of Eq. (10) yields the 
Grad-Shafranov-like equation 



d 2 ip 



dF 



r 1 dip , . 

dr r dr dz 2 ' dtp 



(13) 



with p being determined by equation (10). 

In the limit a <C R the curvature of the stellar surface 
can be neglected on the scale of the accretion region. In 
this approximation U — gz, so that 



p{ip,z) =po(V0 



1 - a 



l/c 



(14) 



where z = corresponds to the top of the crust, 
a = 7 — 1, po is the plasma density at z = 0, h = 
— (dh\p/dz)~ 1 \ z= v = aF/g is the density scale height at 
the footpoints and g is the gravitational acceleration. Note 
that po and h are functions of ip. 

4. ENERGY PRINCIPLE 

We now address the question of the linear stability of 
the polar cap plasma in the approximation described in 
the previous section. We base our analysis on the MHD 
energy principle (Bernstein ct al. 1958). In the present 
section we discuss the energy principle as applicable to 
the polar cap plasmas, generalizing the "intuitive" form of 
the energy principle (Furth et al. 1965; Greene & John- 
son 1968) to include the effects of a strong gravitational 
field. In our discussion we follow the approach of Freidberg 
(1987). 

We start by considering the linearized momentum bal- 
ance equation for a static equilibrium 

p0=F(^) = ijxB+ijxB-Vp-pVC/ (15) 

where £ is the plasma displacement and the gravitational 
potential is unperturbed. Here and in the following the 
tilde denotes perturbed quantities while symbols without 
it stand for the equilibrium quantities. Perturbed current 
J, magnetic field B, density p, and pressure p are respec- 
tively given by Ampere's law, the induction equation, the 
continuity equation, and the equation of state: 



P = 



J = —V x B 

47T 

B = V x (| x B) 
p = -V • (p£) 
— 7pV • £ — £ • Vp. 



(16) 

(17) 
(18) 
(19) 
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We shall assume that the plasma boundary (the crust) is 
rigid and perfectly conducting, so that £ and B vanish 
there. 

The second-order potential energy is given by (Bernstein 
et al. 1958) 



5W = -lj drt-F(t) 



-J x B + - J x B 

c c 



V(tpV • € + € • Vp) + V • (p£)VU 



(20) 



where dr denotes integration over the plasma volume. We 
now transform the above expression for the potential en- 
ergy into a form that allows a physical interpretation. This 
is analogous to what Furth et al. (1965) and Greene & 
Johnson (1968) have done in the absence of gravitational 
field. In order to derive it, we follow the standard approach 
[e.g., Freidberg (1987)], with appropriate modifications for 
the presence of a strong gravitational field. 

Integrating by parts and dropping surface terms (which 
vanish because of the assumed boundary conditions), and 
exploiting the equilibrium relation given by Eq. (10), one 
obtains, after some straightforward algebra, 



4tt 



5W = | / dT { 7 p(V -i-lKg- if + ep(£ • V hip) 2 + 

i ■ [ij x B + 2K g {pi ■ VF) + V(pi ■ VF)] } . (21) 

where K g = (p/-fp)\7U/2. The the second term on the 
RHS, proportional to the Schwarzschild discriminant, can 
be neglected for highly degenerate plasmas in strong mag- 
netic fields because of the smallness of parameter e = 
dlnpp^ 1 /"* /d\np. We shall omit this term in the subse- 
quent analysis (see, however, discussion of finite e effects 
in §9). 

The last two terms in the square bracket on the RHS 
of Eq. (21) can be expanded and recombined to yield an 
alternative form of the energy integral 

5W= l -j dr J|! +7 p(V-£-2K 3 -0 2 



1 ~ p 2 

-J x B + — (£ • VF)/jVF + pV(£ ■ VF) 

c 7p 



Since B • VF = and 



-B J x B = pB • V(£ • VF), 



.(22) 



(23) 



it is clear that the parallel (to B) component of the ex- 
pression in the square bracket in Eq. (22), and hence also 
in Eq. (21), vanishes. Thus the energy integral can be 
expressed in the form 

SW=\jdr{% +lp{V-i-2n g -if 
-€ ± • [iJ x B + 2{p£ ± ■ VF)k 9 + V(p£ ± • VF)] ^24) 

Next, splitting up B into transverse and parallel compo- 
nents, B = Bi + F |b (b = B/F), wc find 



^■JxB = J||£ ± -bxB + F||£ ± -Jxb (25) 
F,| = -F(V • £ ± + 2k c ■ €J + ' VF, (26) 

where k c = b • Vb is the field line curvature. With these 
expressions the potential energy can be expressed in the 
following, "intuitive" form 



6W 



+ 7 p(V ■ £ - 2 Kg ■ O 2 - ^||€± x b B 

-2k • £ ±P {t ± ■ VF)} . (27) 

where k = k c + n g . The first term is the energy cor- 
responding to field line bending, the second is the field 
compression energy, the third is the energy of plasma com- 
pression and buoyancy, the fourth one is responsible for the 
kink instability, and the last one drives the curvature and 
gravitational instabilities. The above expression is a gen- 
eralization of the result of Furth et al. (1965) and Greene 
& Johnson (1968) to the case of a degenerate plasma in 
a strong gravitational field. The above expression reduces 
to these earlier results if the density scale height is much 
longer than the plasma length, i.e., if U <C c 2 s . 

Thus far we have not taken into account the geometry of 
polar caps and the derived energy principle (eq. [27]) ap- 
plies to an arbitrary configuration of plasmas satisfying the 
equation of state (7). We shall now perform a reduction 
of the energy principle as applicable to the configuration 
of our interest. 

First we note that the parallel displacement £|| appears 
only in the nonnegative compressive term. Therefore this 
term can be minimized, in a standard manner, by consid- 
ering the perturbation of the compressive energy due to 
a perturbation of £||. The minimum is found for the dis- 
placement satisfying the Euler-Lagrange equation which 
can be expressed in the form 



B- V 



(V ■ € - 2« fl ■ i) 



= 



(28) 



Integrating along the field line and imposing the line- 
tying boundary condition (£|| [0] = 0) yields 



= 4/' 

p-< Jo 



dl 1 
-B 1 



1 A + V-^-2k 9 -C ± (29) 



where I is the distance along the field line. The integra- 
tion constant A must be determined from the integrability 
condition. Note that p, p — > as z — > z = h/(j — 1) (cf. 
equations [14] and [8]). Therefore, for the kinetic energy 
of parallel motion to be finite, (i.e., J drp^ 2 < oo), this 

requires that 



Jo" gg! (V^i-^^i) 
J? f 



(30) 



where lo = l(z ). This implies, in particular, that £\\(zo) 
0. Since 
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(31) 



and A ^ 0, the compressional energy cannot be minimized 
to zero. This is analogous to that which occurs in flux 
tubes line-tied at both ends (cf. Hameiri et al 1991) and 
is, in fact, caused by the effective line-tying in the low 
density region. 

Setting J|| =0, as appropriate for an untwisted axisym- 
metric field, the potential energy becomes 



SW 



= 3/* ft + £< v ^ + a *-«J J 



+7P" 



Jo" fg_ 

- 2p(€ ± ■ VF)k ■ € ± 1 . (32) 



5. BALLOONING INSTABILITY 

We now limit our attention to the short wavelength (in 
the direction perpendicular to the magnetic field) modes. 
Following the method of Freidberg (1987), we perform the 
reduction of the intuitive form of the energy principle (eq. 
[32]). 

We represent the displacement in the eikonal approxi- 
mation 



SW 



-2p(€ ± -VF)(#6-€x)|. (38) 
From equation (36) it follows that 



B ± = e lb (VX x b) ± + c.c. = (b • VI)b x k ± + c.c. (39) 

Consequently, one finds from equation (38), exploiting the 
hermiticity of the operator F (Bernstein et al. 1958), that 



SW 



= J dr {fci|b • \7X\ 2 + 4(k x b • k^) 2 |X| 



-^(k x b • k ± )(VF x b • kx)|X| 2 | .(40) 

For an axisymmetric system, S = m9 + S(ip,l) where I 
is the distance along the field line. In the absence of a 
toroidal field Bg,B-VS = implies that S = S(ip). Thus 



kx = k e + k n h 



(41) 



where kg — m/r and n = V^>/|VV>| = Vtfj/rB is the unit 
vector normal to the flux surface. Noting that 



£_L = V±e lS + c.c. (33) 

where c.c. stands for complex conjugate; ij ± is assumed to 
be slowly varying while the eikonal S is assumed to vary 
rapidly, i.e., 



kx = VS 



(34) 



is so large that k±a^> 1, where a is the equilibrium length 
scale. 

Note that B • VS = 0. Consequently 



Bx = e lS [W x (77^ x B)]x + c.c. 



(35) 



In the limit of k± — > 00, the dominant terms in SW 
(Eq.[32]) arc those proportional to (V • £ ± ) 2 w |kx ■ f7x| 2 
(second and third term on the RHS of eq. [32]). Therefore 
the minimization of the potential energy, to lowest order 
in l/k±a, requires that k± ■ rj ± = 0, which implies that 
for the zeroth-order displacement 



|b x kx. 



(36) 



where A" is a scalar function. To the next order, the first- 
order displacement can be chosen so that the sum of two 
positive terms involving V • (second and third term on 
the RHS of eq. [32]) is minimized. Since we are interested 
in the limit of (i = 8np/B 2 3> 1, we can instead mini- 
mize only the second of these terms. A simple choice of a 
minimizing displacement is 



2/v n 



£x- (37) 
For this choice, £|| = as well as the third term on the 
RHS of eq. (32) vanish so that 



b x kx = — kgh + k n 8 
.dF 



b x kx • VF = -kgrB 



(42) 
(43) 



the potential energy can be put in form 



SU =4^ / dr{(k 2 e +kl) 



dX 



dl 



+4k 2 n n ( Kn ^f^)\X\ 2 



(44) 



(F' — dF/d^p). Minimization with respect to k n is 
achieved by setting k n = 0. 

The integrand in Eq. (44) involves derivatives of X only 
along the magnetic field but not derivatives of ip. Hence 
ip can be treated as a parameter and the potential energy 
can be expressed in the form 



SW 



J dipSWiip), 



(45) 



where 



dX 



dl 



+4 Kn ( Kn -^l)\X\ 2 



(46) 



Equation (45) implies that SW{ip) < for some tjj is both 
a necessary and sufficient condition for instability (Bern- 
stein et al. 1958). 
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Let us first address the question of stability of a plasma 
confined by a stellar dipole magnetic field. We approx- 
imate F' by AF/ip a where ip a = a 2 B /2 and AF = 
Apo/po, with Apo and po being, respectively, the over- 
pressure and density at the top of the crust. Since n gn = 
(dU/dn)/2c 2 <~ —b r /2h, n cn <~ b r /R, b r ~ r/R (for a 
dipole field) and the density scale height h = c 2 /g is much 
smaller than the neutron star radius R, the second term 
on the RHS of Eq. (46) is destabilizing, if A/3 > a 2 /hR. 

For the interchange modes dX /dl — 0. Because of line 
tying, however, X(0) = with I = being chosen to cor- 
respond to the top of the crust. Expanding into Taylor 
series, we can approximate X ~ I for I < h. From the 
structure of the integral in Eq. (46) we expect that the 
most unstable displacement vanishes for I » h. 

We shall now estimate the integral in Eq. (46). To low- 
est order in a/R, where a is the radius of the polar cap, 
U = c 2 s l/h. We then find that 5W(ip a ) becomes negative 
only if 



R 
h 



(47) 



where A/3 = 8tt Apo / B 2 . The above equation implies that 
for a low overpressure (A/3 < 1), the plasma confined by 
the dipole magnetic field of a neutron star is stable with 
respect to ballooning modes. Equation (47) indicates that 
the plasma may become unstable at high A/3 ^ 1. For 
such a high plasma overpressure, however, the equilibrium 
magnetic field is distorted away from the dipole field and 
has to be computed self-consistently. We consider this case 
in the next section. 

6. A SIMPLE HIGH-A/3 EQUILIBRIUM 

In the limit a <C R the dipole magnetic field in the po- 
lar cap is well approximated by a uniform field. Let us 
therefore assume that prior to the loading of the accretion 
column with plasma, the magnetic field B = Bqz. We 
shall assume that at z — the flux at the is unaffected 
by the accretion, i.e., ip(r,0) = r 2 Bo/2. The loading is 
assumed to occur on flux surfaces with ip < ip a = a 2 B /2. 

Until now, the overpressure profile at the base of the ac- 
cretion region, which determines the function F, has not 
been specified. As a simple example we shall consider the 
following form of F, 



F(V>) = F ext + AF 



1 - 



l+e 



(48) 



for ip < ip a = a 2 B /2 7 with s -C 1; F is a constant F ext 
for tp > ip a . The gradient of F is approximately linear for 
tp < ip a (l — s) and goes smoothly to zero at ip — ip a . In 
the following we shall restrict our attention to the case of 
an overpressure Ap that is small compared to the total 
pressure po, so that AF = Ap /po- In the same approxi- 
mation, the scale height h is, to lowest order, independent 
of ip. 

For such a pressure profile, the Grad-Shafranov equa- 
tion (Eq. [11]) inside the accretion column [ip < ip a (l — e)] 
becomes approximately linear, 

d (ldip\ d 2 ip 4A/3 



We anticipate that because A/3 » 1 and h <SC a, the 
radial excursion of the flux surfaces Ar — J dlb r (< h) 
will be small compared to r. We can therefore neglect the 
r-dependence on the RHS of Eq. (49) and treat x = r / a 
as a parameter. With these approximations, we can then 
solve Eq. (49) by separation of variables, 



1p(r, Z) = ^PoX 2 f(0\ X =r/aX=z/h, 

where / satisfies the equation 



11 

d( 2 



AV(l-aC)-/ 



(50) 



(51) 



(A = 2h\J A/3 / a) . We solve this equation subject to the 
boundary conditions /(0) = 1 and /'(a -1 ) = 0. The first 
of these conditions expresses normalization; the second fol- 
lows from the requirement that the field be vertical in the 
region of vanishing overpressure. By changing variables, 
we can transform the above equation into a modified Bessel 
equation which has the solution 

f(0 = u v (C){AI-4u(0]+Bl4u(0}}, (52) 

where u(() = 2X X (l - a()^ /{I + 2a), v = a/(l + 2a) = 
(7 — 1)/(27— 1), l± v are the modified Bessel functions, and 
A and B are the integration constants. From the condition 
/'(a -1 ) = it follows that B — 0; the boundary condition 
/(0) = 1 then implies that 



/(0 = (i-O- 



2A X 



1+2q 



( 2A X 
~ v \ 1+2q 



(53) 



An example of flux surfaces given by equations (50) and 
(53) is shown in Figure 2. We see that even for relatively 
large A, the radial excursion Ar of flux surfaces is small 
compared to r, in accord with our assumptions, 
z 

0.25 
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0.15 
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0.05 




0.2 



0.4 
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Fig. 2. — Poloidal cross-sections of flux surfaces for A = 1 (solid) 
and A = 0.5 (dashed), with a/h = 10 (r and z are in units of a). 

Let us now consider the case of a/h < A/3 <C a 2 /h 2 
which, as we shall see, is valid at the onset of instability. 
In this case, A <§C 1 and we can expand the Bessel functions 
into a power series. In this approximation we find that 



1 - A' 



r 2 1 - (1 -az/h) 
a 2 (1 + a) (1 + 2 a) 



• (54) 
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With the aid of the above expression we can now demon- 
strate the self-consistency of our earlier assumptions. We 
find that the first term on the LHS of Eq. (49) is smaller 
than the second one as long as p/po > (2ft/a) 2 . Thus for 
h <C a our solution is indeed self-consistent in the region of 
plasma where the most unstable displacement is expected 
to be localized (see §7). 

From (54) the magnetic field components are 



B 

B z :=y B 



(55) 
(56) 



where p = 2A/3hx 3 /-fa. From the above it follows that 
the angle <j> between the magnetic and gravitational fields 
is given by 



tan< 



■ p 



1 



(57) 



From equation (57) one sees explicitly that the field line 
distortion becomes large (tan</> > 1) when p > 1, i.e., 
A/3 > a/ft, as originally observed by Hameury et al. 
(1983). This is, in fact, a simple consequence of balancing 
overpressure Ap/a with magnetic tension kB 2 /4it since 
nh ~ 1 when the field distortion becomes large. 

7. CRITERION FOR THE ONSET OF INSTABILITY 

We now use the equilibrium found in the previous sec- 
tion to evaluate the energy integral. Let us first change 
the integration variable in Eq. (46) from I to z, such that 

dl = dz/b z . For the nearly parabolic pressure profile con- 
sidered earlier, 



+4k„ K n + 




2ft B p a 



(58) 



In the approximation of Eq. (54), we find that to lowest 
order in A 



-(1 + 27&2) 



Po 



b r 
2ft' 



(59) 



Changing the integration variable to s = tan</>, as given 
by equation (57), we then find, with the aid of (55), (56), 
(57) and (59), 



6W(tp) 



p 



1-1/7, 



8Trr 2 -fhB Jo 



,1/7 



+ S 2 



dX 



-(7-1) 



ds 



( S 2 - M 2 )(,s 2 + 2 7 +l) 2 



(60) 



(s 2 + l) 2 

where po — y^y + 1)7(7 — !)■ We observe that the in- 
tegrand in equation (60) is positive for s > pq. Conse- 
quently, the most unstable displacement is expected to 
vanish for s > p a which implies the instability is localized 
to the region where 



P>[^ 
P 



Po- 



(61) 



Since the line- tying in the crust requires that X(s 
p) = 0, we therefore choose as a test function 

Po 



X(s) =X sm[n 



P - Po 

For this test function and 7 = 4/3, as appropriate for a 
plasma in the vicinity of the crust, we find that SW(ip) < 
for p > 11.67 for which it follows that 



A/3 > 7.8 



h' 



(62) 



is a sufficient condition for instability. For p = 11.67, 
equation (61) implies that the instability is localized to 
within 0.93 scale height from the crust. 

The instability criterion (62) is obtained for the equilib- 
rium given by equation (54), derived in the approximation 
of A/3 <C (a/ft) 2 . Thus self-consistency implies that the 
validity of (62) is limited to h/a <C 0.13. This condition is 
satisfied, albeit marginally, for the parameters discussed in 
§2. The instability threshold for more general equilibria, 
with higher values of ft/a, is at present unknown. 

8. MASS OF CONFINED ACCRETED MATTER 

Assuming that the crust is perfectly rigid, the amount 
of mass needed to drive the instability is, from inequal- 
ity (62), 

2 B 2 A{3 

AM = ira 2 - 

8^ g 

= 0.735V-. 

P 

Substituting Eq. (2) in Eq. (63), we have 

2 



(63) 



AM = 3.8 x 10" 13 M P 







B 



10 12 G 



VlO 5 cm J 



Z\ 2/3 /A 



2(> 



56 



5 x 10 8 K 



At a given accretion rate M, the time necessary to accu- 
mulate this mass is AM/M = 33 hr(lO" lo M yr^/M). 
At present we cannot say what, if any, are the observable 
effects of the instability, but the timescale to accumulate 
a critical mass AM is well-suited for monitoring. 

9. DISCUSSION 

Using the MHD energy principle, we have shown that 
the accreted material on the polar cap of a strongly mag- 
netized neutron star is subject to the ballooning instability 
when the overpressure in the vicinity of the crust exceeds 
the magnetic pressure by the factor 7.8a/ft, where a is the 
horizontal pressure gradient length scale (identified with 
the polar cap radius) and ft is the density scale height. 
For typical polar cap parameters, the critical overpres- 
sure is <~ 10~ 4 of the total plasma pressure. The insta- 
bility is localized to within a density scale height of the 
field line footpoint (at the top of the crust) which implies 
that plasma is perturbed at densities (p/po) > 0.33; for 
Po found from equation (1) for a model atmosphere (cf. 
Brown & Bildsten 1999), with a temperature of 7.6x 10 8 K, 
assuming a locally Eddington accretion rate. This im- 
plies p > 10 10 g cm~ 3 . This density is higher than where 
the accreted hydrogen and helium ignite which suggests 
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that the ballooning instability near the threshold for on- 
set will not affect the nuclear burning. It is not clear, 
however, whether the instability can degrade the confine- 
ment enough to prevent a further accumulation of mass 
well above the instability threshold. 

In our stability analysis we have omitted the contri- 
bution of the Schwarzschild term to the MHD potential 
energy arguing that it is small in highly degenerate plas- 
mas of uniform composition immersed in strong magnetic 
fields. To test the validity of this approximation we have 
performed a calculation along the same lines as in §4 and 
§5, including the buoyancy contribution, for e = 4 x 10~ 4 , 
for a realistic atmosphere (cf. Brown & Bildsten 1999). 
We found that this inclusion increases the numerical fac- 
tor on the RHS of inequality (62) by 17% for B = 10 12 G 
and by less than 1% for B > 5 x 10 12 . Thus the effect of 
buoyancy is indeed weak for magnetic fields of interest, as 
long as e is as small as assumed. 

The latter assumption, as well as that of composition 
uniformity, is violated in the electron-capture layers. The 
location of these layers is not precisely known. Haensel & 
Zdunik (1990) found that electron capture layers closest 
to the top of the crust occur at densities 10 10 g cm~ 3 and 
8 x 10 10 g cm~ 3 . Since the instability found in the present 
paper is localized within a density scale height from the 
crust which for the model atmosphere implies localization 
between densities 10 10 and 3 x 10 10 g cm~ 3 , we expect that 
the electron capture effects would not significantly modify 
our conclusions. 

In our analysis the crust was modelled as a rigid bound- 
ary in which the magnetic field is anchored. For this ap- 
proximation to be justified the lateral distortion of the 
crust in response to the overpressure must be much smaller 
than the magnetic field line distortion in the overlying 
fluid. This is true when A/3/ j3 exceeds a s Aa/a where 



Aa <~ (ft. 2 /a) A/3 is the displacement of flux surfaces in 
the ocean (cf. §6), i.e., when (3 > (a/h) 2 /a s ; here a s p is 
the crust shear modulus. Estimates (Pandharipande et al. 
1976; Strohmayeret al. 1991) find a s ~ l(T 3 -l(r 2 , which 
for the characteristic density scale height and polar cap ra- 
dius (cf. §1) implies that (3 > 10 4 - 10 5 . For the model 
atmosphere, this condition is satisfied for B < 5 x 10 12 G. 

The imperfect rigidity of the crust may affect the esti- 
mate of the accreted mass necessary to destabilize the bal- 
looning mode (eq. [63]). Ushomirsky et al. (2000) found 
that the crust tended to sink in response to a lateral den- 
sity discontinuity. While their findings are not directly 
applicable to our case, it is nevertheless conceivable that 
the crust could deform sufficiently to significantly reduce 
the applied overpressure. The mass required to drive the 
instability would then be higher than our estimate. 

The effect of the instability is an enhanced cross-field 
transport of matter and a resulting decrease in the plasma 
confinement time in the polar cap. To find the magnitude 
of this anomalous transport requires analysis of nonlinear 
effects which are beyond the scope of this paper. The anal- 
ysis of these effects and, in particular, the question of the 
existence of a steady state, in which the transport is suffi- 
ciently rapid to balance the influx of the accreted matter, 
are relegated to future publications. 
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